Simulations of liquid crystal hydrodynamics 
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We present a lattice Boltzmann algorithm for liquid crystal hydrodynamics. The coupling be- 
tween the tensor order parameter and the flow is treated consistently allowing investigation of a 
wide range of non-Newtonian flow behavior. We present results for the effect of hydrodynamics on 
defect coalescence; of the appearance of the log-rolling and kayaking states in Poiseuille flow; and 
for banding and coexistence of isotropic and nematic phases under shear. 
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There is growing interest in obtaining a fundamental 
physical understanding of the flow properties of liquid 
crystals, polymer melts, and droplet suspensions. The 
hydrodynamics of such complex fluids can be complicated 
and very different from that of simple liquids because of 
the coupling between the microscopic structure and the 
velocity fields imposed by the flow ^ . Examples include 
shear-thinning and thickening and non-equilibrium phase 
transitions such as banding under shear |^-^]. 

To understand the physics underlying such flow prop- 
erties it is helpful to develop simulation techniques to 
probe the hydrodynamics of complex fluids. This has 
proved difficult because of the diverse length and time 
scales involved. Microscopic approaches, such as molec- 
ular dynamics, provide the most faithful representation 
of the microscopic physics. However they are not usually 
able to probe hydrodynamic time scales. 

One attempt to surmount this problem has been the 
development of lattice Boltzmann simulations Q . These 
solve the hydrodynamic equations of motion while in- 
putting sufficient, albeit generic, molecular information 
to model the important physics of a given fluid. This is 
often done by imposing a Landau free energy functional, 
so the fluid evolves to a known thermodynamic equilib- 
rium 

Most of the results thus far using lattice Boltzmann ap- 
proaches have concentrated on the properties of binary 
fluid mixtures. Considerable progress has been made in 
understanding the effect of hydrodynamics on domain 
growth in 2- and 3-dimensions Q . Investigations of flow 
have been more limited although there has been some 
work on binary solutions under oscillatory shear, on the 
flow of binary mixtures in porous media, and on am- 
phiphilic fluids ||^. 

In this paper we describe a lattice Boltzmann algo- 
rithm for liquid crystal hydrodynamics. The aim is to 
investigate the wide range of physical phenomena which 
result because the director field couples to the flow. First 
we assess the effect of hydrodynamics on defect coales- 
cence, a process important for phase ordering. This illus- 
trates how the director configuration induces flow. We 
then examine the orientation of the director in Poiseuille 
flow. In addition to a steady state director configuration 



for slower fluid velocities, the director can exhibit tum- 
bling and chaotic orbits in rapid flows. Finally, we study 
the system under shear and demonstrate how the liquid 
crystal undergoes a non-equilibrium phase transition to 
a banded state where the bands are different phases co- 
existing at different strains but a unique stress [||j3). 

We emphasize that almost all previous work has ei- 
ther ignored the flow, imposed a flow and worked out 
the director configuration ignoring back-flow effects, or 
coupled states of the director with states of flow using 
arguments based on the stability of interfaces. Here re- 
sults are obtained self-consistently, by simulating the full 
hydrodynamics equations coupled to a description of the 
liquid crystal based on a tensor order parameter. 

First we present the relevant equations of motion and 
the extensions of the lattice Boltzmann approach needed 
to solve them. A major difference from more simple 
fluids is that liquid crystals are described by a ten- 
sor order parameter Q (related to the director n by 
Qai3 = {natip — ^5ap))- Their equilibrium properties 
can be described by the Landau-de Gennes free energy 
functional 
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where (j) is the liquid crystal concentration and T the 
temperature. The exact form of the coefficients of the 
bulk free energy terms is not important and for simplic- 
ity we write them in terms of a single parameter 7 pO[ . 
Similarly we restrict ourselves in the first instance to a 
single elastic constant k. 

The order parameter is not conserved. It evolves ac- 
cording to the convection-diffusion equation P,p|,pd| 

{dt + 9„Ua)Q = S(W, Q) + {T/kBT^mq) (2) 

where u is the bulk fluid velocity and F is a diffusion con- 
stant. Of the terms on the right-hand side of Eqn. (^ 

S(W, Q) = (AD + n)(Q + ll) + (Q + ^I)(AD - fl) 



-2A(Q + il)Tr(QW) 



(3) 
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accounts for the rotation of the nematic order parameter 
driven by the symmetric D, and antisymmetric fl, parts 
of the velocity gradients Wap = dpUa- The parameter 
A is related to the aspect ratio of the polymer molecule 
or alternatively can be viewed as a phenomenological pa- 
rameter to correct for inaccuracies of quadratic closure 
HI. The molecular field, 
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describes the evolution to equilibrium in a way analogous 
to Model A. 

The flow of the liquid crystal fluid obeys the continuity 
and Navier-Stokes equations 



dtp + dapua = 0, 



(5) 



pdtUa + pUpdpUa = dpTap + dpCTap + 

2pTf/'i{di3{5al3 + idn(Jap)dc_UQ + daUp + dpUa) (6) 

where 

Oap = —pTSap — SHap — 

K{daQcsdpQ<:s - dxQc&dxQ^&dap/'i) (7) 
T^p = HQ - QH (8) 

are the symmetric and antisymmetric contributions to 
the non-dissipative part of the stress tensor respectively. 

A lattice Boltzmann scheme which reproduces equa- 
tions (0) to (jsj) to second order can be defined in terms 
of two distribution functions fi (x) and (x) where i la- 
bels lattice directions from site x. Physical variables are 
related to the distribution functions by 



a = E/^e,„, Q = ^G,. (9) 



The distribution functions evolve in a time step At 
according to 

Mx + e,At, t + At)- Mx, t) = ^ [C/,(x, t, {/,}) + 
Cf,ix + e,At,t + At,{f*})] 

G,(x + e,At, t + At)~ G,(f , = ^ [CG^{x, t, {G,;})+ 



CG^{x + e,At,t + At,{G*})] 



(10) 



where the collision operators are taken to have the form 
of a single relaxation time Boltzmann equation, together 
with a forcing term 

--{h{x, t) - f^(x, t, {/.})) t, {/,}), 

CGi(x,t, {Gi}) = 

- — (G,(x,t)-G°(f,t,{Ga)) + /i,;(x,t,{G,;}). (11) 
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f* and G* in equations (|10|) and ( |llD are first order ap- 
proximations to fi {x + CiAt, t + At) and G^ [x + e,; At, t + 
At). They are introduced to remove lattice viscosity 
terms to second order and they give improved stability. 

The form of the equations of motion and thermody- 
namic equilibrium follow from the choice of the moments 
of the equilibrium distributions /f and G^ and the driv- 
ing terms pi and hi. is constrained by 



i 

^ ^ fi SiaGip 



^ ^ fi ^ia — PUa, 

i 

-0-ap + PUaUp 



(12) 



where the zeroth and first moments are chosen to im- 
pose conservation of mass and momentum. The second 
moment of f'^ controls the symmetric part of the stress 
tensor, whereas the moments of pi 

= 0, ^PiCia = dpTap, ^PiCiaBip = (13) 



impose the antisymmetric part of the stress tensor. 

For the equilibrium of the order parameter distribution 
we choose 



^G? = Q, ^G?e,„ = Qi 

i i 
E GiCiaCip = QUaUp. 



(14) 



This ensures that the fluid minimizes its free energy at 
equilibrium and that it is convected with the flow. Fi- 
nally the evolution of the order parameter is most con- 
veniently modeled by choosing 

J2 = {D/ksTmiQ) + S(W, Q), 

i 

i i 

Conditions (p^-(p^ can be satisfled by taking /f , G°, 
hi, and pi as polynomial expansions in the velocity as is 
usual in lattice Boltzmann schemes Q. A second order 
Chapman-Enskog expansion for the evolution equations 
( p^ ) incorporating the conditions ([l^)-(p^ leads to the 
equations of motion (||), (||), and (j^ 

The alternative Ericksen-Leslie-Parodi (ELP) equa- 
tions of liquid crystal hydrodynamics are written in terms 
of the evolution of the director field n Q rather than the 
tensor order parameter Q. We use the latter approach 
for two reasons. First, the motion of disclination lines 
(points in two dimensions) is explicitly included and it is 
particularly interesting to assess the effect of flow on this 
dynamics. Second, we are interested in examining phase 
transitions induced by shear and so need to describe both 
the isotropic and nematic phases. The ELP formalism 
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describes only the nematic phase at a fixed amphtude of 
order parameter. If we restrict the order parameter to 
be uniaxial with fixed amplitude, the Q equations of mo- 
tion can be reduced to the ELP equations. However, due 
to our choice of coefficients in, for example the free en- 
ergy, we have specialized to three independent viscosity 
coefficients rather than the five in the ELP formalism. 
A lattice Boltzmann method for implementing the ELP 
equations is being studied elsewhere 

If Fig. H we show the flow fields induced when two 
disclinations of opposite sign mutually annihilate, a pro- 
cess important in phase separation. The director field 
induces a fiow in the fluid in the form of vortices ac- 
companying the disclinations Q . Even after the defects 
have annihilated, they leave vortices behind in the fluid 
which then gradually decay away. The vortices induced 
in the fluid speed up the annihilation process, but do not 
change the t^^'^ power law ||l^,|l5| for annihilation. This 
is because the vortex-vortex interaction is of the same 
form as the disclination-disclination interaction. How- 
ever, in a system with many disclinations the interaction 
between annihilating pairs is likely to be much stronger 
when hydrodynamics is present as the fluid vortices left 
by the annihilation process are at a considerable distance 
from the annihilation site. We are currently examining 
the effect of hydrodynamics on the phase ordering. 

A second example of the coupling between the director 
field and the flow is shown in Fig. || where a Poiseuille flow 
is imposed on the liquid crystal with the director field 
pinned perpendicular to the flow direction at the bound- 
aries. The director field couples to the shear component 
of the flow leading to a rich variety of possible dynamical 
states. For slow flows a static director field configuration 
is induced by the fiow (Fig. §(a)). We note that in or- 
der to match the director fields at the midpoint of the 
flow, the director orients parallel to the flow, not per- 
pendicular as has been previously suggested js), as this 
configuration has the lower elastic energy |l6j. For faster, 
more rotational flows, or alternatively smaller A, the "log 
rolling" state, where the director points perpendicular to 
the shear plane, has a lower viscosity and is more stable 
in the bulk. However, as the director is pinned in the 
shear plane at the boundary it must match the bound- 
ary configuration onto the log rolling state in the interior. 
The system is unable to do this with a steady state con- 
figuration and the director tumbles and kayaks (rotates 
in and out of the plane) in the interface region, as show 
in Figure §(b). While the velocity profile remains nearly 
parabolic, as would be expected for Poiseuille flow in a 
simple fluid, the effective viscosity is strongly influenced 
by the state of the director. 

Shear acts to orient the director field with respect 
to the flow and in this sense is analogous to a mag- 
netic fleld in spin systems. As such, it can shift the 
transition from the isotropic (I) to nematic (N) phases 
P,PHr^ . The resulting non-equilibrium phase diagram in 
the shear-stress-(l — 7/8) plane (see Eqn. |l|) consists of 
a line of first order phase transitions connecting to the 



equilibrium transition at zero stress and terminating at 
a non-equilibrium critical point. This is illustrated in 
Figure ||(a) which shows the behavior of Si , the largest 
eigenvalue of Q, as a function of the shear stress. As 7 
is increased, the jump in the order parameter decreases 
and disappears at a non-equilibrium critical point. 

In an experiment, or in our simulations, we do not have 
direct control over the shear stress (i.e. the total off di- 
agonal stress, including elastic and viscous terms), but 
rather control the average strain rate (the relative veloc- 
ity of the walls). The coexistence region is just a line on 
the stress-temperature phase diagram, but it corresponds 
to a finite region on a strain-temperature diagram. In 
fact the two different coexisting phases have a different 
strain rate. This phenomenon is known as shear banding 
and results in a plateau in the stress-strain curve shown 
in Figure ||(b). Once coexistence is reached, increasing 
the strain rate does not increase the stress, it just con- 
verts more of the system from the I to N phase. Once all 
the system is in the N phase, the stress again starts to in- 
crease. Our results are in agreement with those obtained 
in Ref. |^,^. There an interface was imposed on a one- 
dimensional system and its stability was used to infer the 
state of the system. Here, shear banding is spontaneous 
and work is in progress to investigate the way in which 
the bands form. 

In conclusion, we have derived and implemented a lat- 
tice Boltzmann algorithm for liquid crystal hydrodynam- 
ics. This opens the way to investigate a wide range of 
physical phenomena which result from the coupling be- 
tween the director field and the flow, many of which have 
been examined only indirectly or with severe approxima- 
tions in the past. 

We thank G. Gonnella, C. Care and P. Olmsted for 
useful conversations. This work was funded in part by 
grants from the EPSRC. 
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FIG. 2. Two different states in Poiseuille flow, where the 
lines represent the director orientation (eigenvector corre- 
sponding to largest eigenvalue of Q) projected down onto the 
xy-p\ane, and shading represents the amplitude of the order 
parameter (largest eigenvalue). Flow is from top to bottom, 
and the walls are at the left and right. At the walls, the 
director is aligned perpendicular to the boundary, (a) A sta- 
ble conflguration at low flow, (b) Snapshots of an oscillating 
configuration where the central region is in the "log-rolling 
state" (director perpendicular to the plane) and the bound- 
ary region consists of a transition from a configuration in the 
shear plane to a "tumbling" and "kayaking" region (director 
rotating in and out of the plane) interfacing to the central 
log-rolling state. 



FIG. 1. (a) Director field for ±1/2 disclinations moving to- 
wards annihilation, (b) Fluid velocity induced by the discli- 
nations at the same time as the conflguration shown in (a), 
(c) Vortices left in the fluid velocity fleld after the defects have 
annihilated, (d) Separation of the disclinations as a function 
of time for systems of sizes L — 64, 128 and 256 (left to right) 
with (solid) and without (dashed) hydrodynamics. The discli- 
nations start out at rest at a separation of L/2 — 1. Defects 
in systems with hydrodynamics annihilate faster. 
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FIG. 3. (a) Amplitude of the order parameter (largest 
eigenvalue of Q) as a function of the shear stress for 7 = 2.65 
(squares), 2.6 (diamonds), 2.55 (triangles) and 2.5 (stars). 
The equilibrium transition occurs at 7 « 2.7. (b) The 
stress-strain curve shows a discontinuity at the point where 
the shear stress induces a first-order, non-equilibrium phase 
transition from the weakly ordered paranematic I phase to 
the strongly ordered nematic N phase. 
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